Green-Function-Based Monte Carlo Method for Classical Fields Coupled to Fermions 
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Microscopic models of classical degrees of freedom coupled to non-interacting fermions occur in many dif- 
ferent contexts. Prominent examples from solid state physics are descriptions of colossal magnetoresistance 
manganites and diluted magnetic semiconductors, or auxiliary field methods for correlated electron systems. 
Monte Carlo simulations are vital for an understanding of such systems, but notorious for requiring the solution 
of the fermion problem with each change in the classical field configuration. We present an efficient, truncation- 
free O(N) method on the basis of Chebyshev expanded local Green functions, which allows us to simulate 
systems of unprecedented size N. 

PACS numbers: 02.70.Ss, 71.15.-m, 75.47.Lx 



The numerical simulation of quantum lattice models is a 
key tool in solid state research and many other fields of 
physics. One class of problems, which is notoriously dif- 
ficult to study, are fermions coupled to classical degrees of 
freedom. Such microscopic models can arise if parts of a 
complex system are approximated classically. A prominent 
example is the double-exchange model, which describes the 
ferromagnetism of mixed-valence manganites on the basis of 
classical <2g -spins whose orientation affects the kinetic energy 
of e g valence electrons. HI Q. S Another example are Mn- 
doped (III,V) semiconductors, where itinerant holes trigger 
a ferromagnetic ordering of the Mn spins. J3l A completely 
different route that leads to a coupling of fermions and classi- 
cal degrees of freedom are the auxiliary field methods, which 
tackle the problem of interacting fermions. Here, a Hubbard- 
Stratonovich transformation is used to decouple the two-body 
interaction into non-interacting fermions in an auxiliary field, 
which is summed over with Monte Carlo methods. Sla] 

For all these systems the cause of the numerical difficulty is 
the requirement for a solution of the non-interacting fermion 
problem whenever the classical field is varied in a Monte 
Carlo simulation. We propose an efficient local-update algo- 
rithm, which obtains the change in fermionic energy directly 
from a few local Green functions. These Green functions can 
easily be calculated by Chebyshev expansion, superseding es- 
timates of the density of states and thus trace calculations. We 
illustrate the efficiency of the approach with simulations of the 
double-exchange model. 

The models we consider in this work are of the general form 
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where are fermion creation (annihilation) operators at 
lattice site i, and is a classical field with one or more 
components at each site. For example, the double-exchange 
model [7, 8] is given by the Hamiltonian 



H — tij c i c j > 



(2) 



where the summation is over nearest-neighbor sites and the 
hopping tij depends on the orientation {0i,<fii} of classical 



local spins at each site, 

£,„• = cos n 1 cos 



Bi+Oj ■ fa—a 

i cos — y 2 - sin — y- 



(3) 



This complex matrix element is one for ferromagnetically 
aligned spins and vanishes for anti-ferromagnetic alignment. 
At low temperature the system favors ferromagnetism, since 
it can gain kinetic energy. 

The thermodynamics is described by the partition function 



Z = Tr c Tr f exp[-/3(i/(0) - nN)] 



(4) 



and its derivatives. Here the traces Tr c and Trf sum over the 
classical and fermionic degrees of freedom, respectively. The 
fermionic trace can be rewritten in terms of the single-particle 
eigenvalues ej of H, 

Z = Tr c cxp[-S eff (0)], (5) 

such that the grand potential of the fermions times j3, 

S eS ($) = -J2 lo S i 1 + exp[-P{u($) - M)]) , (6) 

i 

defines an effective Euclidean action for the classical degrees 
of freedom. The second trace over the classical field can then 
be calculated with a standard Monte Carlo sampling, where 
the weight of a configuration <j> is given by 



P(<t>)=en>[S QS {4>)]lZ. 



(7) 



However, a non-trivial problem remains: To calculate P{4>) 
we need to know the spectrum {e;} of the non-interacting 
fermion system, or at least its change under a proposed Monte 
Carlo update <j> — * <f>' . 

There are different solutions to this problem: We could use 
brute force and calculate all eigenvalues of (tp) whenever p 
is modified. For local update schemes this is, of course, very 
expensive and imposes severe restrictions on the accessible 
system sizes. As a way out, hybrid approaches Jgl, T/J, 11, 121 
have been suggested, where updates of the whole field config- 
uration are calculated with an approximate dynamics. Then, 
the solution of the full fermion problem in the acceptance step 



is required less frequently. However, if the approximate dy- 
namics does not closely match the exact one, the acceptance 
rate drops markedly, in particular for increased system size. 
Therefore, these approaches crucially depend on the quality 
of the approximate action. 

Staying with local updates of the classical field one can try 
to optimize the calculation of the fermion spectrum. Motome 
and Furukawa 

EI El El 

suggested a Chebyshev expansion 
of the fermion density of states, which can be calculated with 
an effort proportional to the square of the system size TV. A 
further modification [|16ill7ll . involving several truncations in 
the moment calculation, reduced the effort to order TV. 

In another recent approach iflitl the evolution of the eigen- 
values of Aij ((/>) under small local changes of is tracked us- 



ing special techniques for low -rank matrix updates [19]. Then 
again, the full solution of the fermion problem is required only 
occasionally. In the best case this leads to TV log TV scaling. 

In the present work we combine ideas from the last two 
approaches and directly calculate the change of the fermion 
density of states using a few real space Green functions. Re- 
lying on Chebyshev expansion these can be calculated with an 
effort proportional to the system size TV. Without any trunca- 
tions we arrive at an order- TV algorithm. 

Let us start from the Hamiltonian H with the Hermitian 
hopping matrix A and ask how the spectrum changes, when a 
local modification A is added. Given a(E) — A — Et and its 
inverse G(E), i.e. the Green function with G(E)a(E) = 1, 
the spectrum of A + A follows from 



(A 
HE) 



AM 
AM 



0. 



(8) 



Right multiplication with G(E) yields 

G(E)(a(E) + AM = (1 + G(E)A)\ip) = . (9) 
A is given by those values of E where 



The spectrum of A 
the determinant 



d(E) := det(l + G(E)A) 



(10) 



vanishes. Recalling introductory lectures on Green functions 
(see e.g. Ref. I20I) we note that this TV-dimensional determi- 
nant reduces to one with a dimension equal to the rank of A. 
For local Monte Carlo updates this is a small number. If we 
change, for instance, the on-site potential, A has a single non- 
zero matrix element and we merely need the local Green func- 
tion of the corresponding site. For the double exchange model 
a single spin flip affects the hopping between the site and its 
nearest neighbors. Independent of the system size TV or the 
space dimension, Eq. ( TTOb then reduces to a 2 x 2 problem, 
i.e., we need only four Green functions connecting the site 
with its environment and both to themselves. 

Before we go into the details of calculating d(E), let us 
further analyze the meaning of this quantity. Going back to 
Eq. (|9]l we have 



d(E) = det[G(E)(a(E) + A)] 
= det[G(E)]det[a(E) + A] 
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FIG. 1: (color online), (a) Labeling of the nearest neighbor sites 
on a cubic lattice, (b) Comparison of Im log(d(_B + i e)) calculated 
from exact eigenvalues e 4 , d i and with Chebyshev expansion of order 
M = 512 (e = 0.0625) for a 6 3 site sample. 



This can be expressed in terms of the eigenvalues {e^} of A 
and {€[} of A + A, 



d ( s )=n^bn«- s ) 



(12) 



An important trick of our new approach consists of going to 
the complex plane, E — > z := E + i e, and taking a logarith- 
mic derivative. We now observe that d(E) determines exactly 
what we need for the Monte Carlo update: the change in the 
density of states going from A to A + A, i.e., from <f> to <j>' , 



Ilmlim^fM 

7T dz 



= — Im lim 



1 



z e — z 
= J2S(e i -E)-S(e' i -E) 

i 

= p{E)-p'{E). 
The change in the effective action then reads 



(13) 



Ses{<t>') - S e «{4>) = 

log(l + e-^ E -^)(p(E) p'(E))dE (14) 

= i + e llE-») Imlimlog(d(g + i £ ))d£;, 

where in the last line partial integration led to an integral over 
the Fermi function. The key role of d(z) has already been 
noted earlier 12 ll 12211 . but the evaluation of Eq. ( TTOb becomes 
feasible only if we can restrict ourselves to a minimal set of 
Green functions and an efficient, direct method for their cal- 
culation. 

Let us now explain this main part of our new approach for 
the double exchange model on a cubic lattice. Here, a local 
update consists of rotating a single spin at a site o. This modi- 
fies the matrix element ty between o and its nearest neighbors 
to the north, east, south, and so on. Labeling the sites accord- 
ing to Fig.fTfa), a naive evaluation of Eq. (TTOb requires all 7 x 7 



3 



Green functions Gij(z) with i,j G {o,n,e, s,w 7 t,b} MM- 
However, we can do much better observing that 



d(z) = det(l + G(z)A) 

= [1+ E A ooG oj (z)][l+ ]T ^ojG 0O {z)] 

-G 00 (z)[ A jo A ok G kj {z)] 



(15) 



can be expressed in terms of only 2x2 Green functions, 

d{z) = [1 + G ov (z)][l + G vo (z)\ - G 00 (z)G vv (z) , (16) 
which connect the original site o and the environment state 



= A\o) = J2 *io\3) 



(17) 



jen.n. 



All four Green functions can be calculated easily with 
the Chebyshev expansion approach outlined in a recent re- 
view 1I24I1 . In a nutshell, diagonal elements Gu{z) are ex- 
panded in terms of the Chebyshev polynomials of first and 
second kind, T m and U m , respectively, 



G lt (E 
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_1 [i m exp[— i m arccos(i?/s)] 



. (18) 

Vs 2 - E' 2 

The expansion coefficients /i m are Chebyshev moments mod- 
ified by appropriate kernel factors, which improve the conver- 
gence of the truncated series and damp Gibbs oscillations, 

sinh[A(l - m/M)} 



Mm = (i\T m (H/s)\i) 



sinh A 



(19) 



The scaling factor s ensures that the spectrum of the Hamilto- 
nian H/s falls within the domain of the Chebyshev polynomi- 
als [— 1 , 1] . For the double exchange model we can choose s to 
be a little larger than the bare bandwidth in the ferromagnetic 
case, s > 6. The kernel parameter A regulates the resolu- 
tion of the method versus the damping of Gibbs oscillations. 
A good value is A = 4. The resolution at which the Green 
function is approximated is given by e = Xs/M. The limit 
e — > thus corresponds to infinite expansion order M. In a 
numerical simulation, of course, M is always finite and we 
need to extrapolate data for different M to obtain the limiting 
value of S e ff(4>') — S f e ff(0)- Since in Eq. ( fT4l we integrate a 
function of Gij (z) over the Fermi function, the maximal res- 
olution should be better than the thermal broadening of the 
Fermi step, which is of the order of 1 / (3. Low temperatures, 
therefore, require higher expansion orders. 

The most time consuming step of the whole simulation is 
the calculation of the moments (i\T m (H / s)\i) . Using the re- 
cursion relation 




FIG. 2: (color online). Main panel: Magnetization versus temper- 
ature for the double-exchange model on 3D clusters with periodic 
boundary conditions (iV = L 3 ). Inset: Binder ratio Ua showing a 
crossing near T ~ 0.134. 



it reduces to sparse matrix-vector multiplications, the cost of 
which scales linearly with the system size N. With a further 
trick based on 



T, 



2m+i 



— 9T T — T 



with i = 0, 1 



(21) 



T m (x) = 2xT m ^x(x) - T m - 2 {x) 



(20) 



we obtain two moments per matrix-vector multiplication. 
Moreover, half of the moments vanish due to the relation 
\v) = A\o) and the special structure of the Chebyshev re- 
cursion (1201) . 

At this point we should also note the advantage over ap- 
proaches based on a full expansion of the density of states: 
For the latter the moments are given by traces, fj, m ~ 
Tr{T m (i//s)}, instead of simple expectation values. Unless 
truncated or otherwise approximated [16], their calculation re- 
quires 0(N 2 ) operations. 

So far we have discussed only the diagonal Green functions 
Gu. In Ref. |24j we showed that symmetric Green functions 
Gij — Gji can be expanded in the same manner. However, 
in the double exchange model the spins induce local magnetic 
fields which break this symmetry. We therefore derive the off- 
diagonal Green functions G ov and G vo from the two diagonal 
functions G + V:0 + v and G +i V . +i v . With all required mo- 
ments at hand we can evaluate the sums in Eq. ( fl~8l with fast 
Fourier methods, calculate d(E) with Eq. ( fl~6l ). and finally in- 
tegrate over the Fermi function to obtain the change of S e s. 
In Fig.QJb) we show a typical example of lm\og(d(E + ie)) 
and compare the expansion with the exact result from a full 
diagonalization. 

Having explained the technical details of the approach let 
us now illustrate its applicability with a few results for the 
double-exchange model (f2]) at half filling, fj, = 0. In the main 
panel of Fig. [2] we show the magnetization as a function of 
temperature, where the latter is measured in units of the max- 
imal hopping amplitude t = 1. As expected, we observe a 
phase transition to a ferromagnetically ordered phase below 
T w 0.14. A closer inspection based on the Binder parame- 
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FIG. 3: (color online). Structure of the 3D spin field at two different 
temperatures visualized with stream lines. 



terQ 



Ua = 1 



3(m 2 



(22) [io] 



yields the estimate T c sa 0.134, see the inset of Fig. [2] This 
agrees with previous estimates llll[|l5U26ll of the critical tem- 
perature, which range between 0. 128 and 0. 139. 

The data in Fig.|2]is based on expansions of order M = 256 
and averages over 6000 to 20000 Monte Carlo steps in the crit- 
ical region, where one step corresponds to N spin flips. It is 
the low resource consumption which allows for these far more 
precise calculations compared to previous studies 11151 12611 . 
which were based on Chebyshev expansions of order M s=s 
20. Note also that we can handle much larger systems with 
N = 30 3 sites, i.e. a complex matrix dimension of 27000. 

To visualize the spatial structure of the spin field we imag- 
ine it as a liquid flow and draw curves tangent to the velocity 
field. These stream lines are regular and parallel in the or- 
dered phase, but quite irregular and swirled near and above 
the phase transition, see Fig. [3] Of course, this visualization 
method fails for truly disordered spin fields. 

In summary, we have presented an efficient local update 
scheme for Monte Carlo simulations of classical fields cou- 
pled to fermions. At the core of the approach is an expres- 
sion which relates the change of the fermion spectrum to a 
few local Green functions. These can be calculated easily 
with Chebyshev expansion. Compared to similar expansion 
approaches a full trace over the fermion system is spared, 
which directly leads to a fast and precise order- ./V algorithm. 
Possibly our method can be further accelerated using the ap- 
proximations inherent in the truncated polynomial expansion 
method HHO- The calculations presented were performed 
on the TeraFLOPS cluster of the Institute for Physics at Greif- 
swald university. 
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